/*
 * case_1D1t_GA.cpp
 *
 *  Created on: 11-09-2012
 *      Author: mith
 */

#include "case_1D1t_GA.h"
// MY HEADERS
#include"mros-v42.h"

namespace case_1D1t_GA
{

global glob;

// *******************************************************************************************  //
// *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY *** START OF THE COPY ***
// *******************************************************************************************  //

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double sdelta;
	double J;
	double t;
	double U;
	unsigned int L;
};

struct variables {
	double BDelta;
	double chi;
	double m;
	double d;
	double mu;
	double Energy;
};

/*
 *********************************************************
 ***************** Small useful functions ****************
 *********************************************************
*/

void clear(parameters &p) {
	p.sdelta = p.J = p.t = p.U = p.L = 0;
}

void clear(variables &v) {
	v.BDelta = v.chi = v.m = v.d = v.mu = v.Energy = 0;
}

void copyTo(variables &v, variables &q, double d=1.) {
	q.BDelta = v.BDelta*d;
	q.chi = v.chi*d;
	q.m = v.m*d;
	q.d = v.d*d;
	q.mu = v.mu*d;
	q.Energy = v.Energy*d;
}

void addTo(variables &v, variables &q, double d=1.) {
        q.BDelta += v.BDelta*d;
        q.chi += v.chi*d;
        q.m += v.m*d;
        q.d += v.d*d;
        q.mu += v.mu*d;
        q.Energy += v.Energy*d;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters
	variables v; // variables
	double n, r, w; // contractions

	void refresh() {
		n = 1.-p.sdelta;
		r = n/2. + v.m;
		w = n/2. - v.m;
	}

public:
	Func(const parameters &p, const variables &v) : p(p), v(v) { refresh(); }

	void setPar(const parameters &p) { this->p = p; refresh(); }
	void setVar(const variables &v) { this->v = v; refresh(); }

	double gt();
	double gs();
	double gBDelta();
	double gm();
	double dgtdm();
	double dgtdd();
	double dgsdm();
	double dgsdd();
	double mygamma(double kx, double ky);
	double eta(double kx, double ky);
	double epsilon(double kx, double ky);
	double BDelta_d();
	double BDelta_af();
};

double Func::gt() {
	double part1 = (n-2.*v.d)/(n-2.*r*w);
	double part2 = sqrt((1.-w)*(1.-n+v.d)/(1.-r)) + sqrt(w*v.d/r);
	double part3 = sqrt((1.-r)*(1.-n+v.d)/(1.-w)) + sqrt(r*v.d/w);
	return part1*part2*part3;
}

double Func::gs() {
	return pow((n-2.*v.d)/(n-2.*w*r),2.);
}

double Func::gBDelta() {
	double part1 = (n-2.*v.d)/(2.*(n-2.*r*w));
	double part2 = sqrt((1.-w)*(1.-n+v.d)/(1.-r)) + sqrt(w*v.d/r);
	double part3 = sqrt((1.-r)*(1.-n+v.d)/(1.-w)) + sqrt(r*v.d/w);
	return part1*(part2*part2 + part3*part3);
}

double Func::gm() {
   return (n-2.*v.d)/(n-2.*w*r);
}

double Func::dgsdm() {
	return -8*v.m*gs()/(n-2*r*w);
}

double Func::dgsdd() {
	return -4*gs()/(n-2*v.d);
}

// Using Mathematica
double Func::dgtdm() {
	return -(((sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.)) +
	          sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.)))*(-2*(-v.m + n/2.) + 2*(v.m + n/2.))*
	        (sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.)) +
	          sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.)))*(-2*v.d + n))/
	      pow(-2*(-v.m + n/2.)*(v.m + n/2.) + n,2)) +
	   ((sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.)) +
	        sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.)))*
	      ((-(((1 + v.d - n)*(1 - v.m - n/2.))/pow(1 + v.m - n/2.,2)) -
	           (1 + v.d - n)/(1 + v.m - n/2.))/
	         (2.*sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.))) +
	        (v.d/(-v.m + n/2.) + (v.d*(v.m + n/2.))/pow(-v.m + n/2.,2))/
	         (2.*sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.))))*(-2*v.d + n))/
	    (-2*(-v.m + n/2.)*(v.m + n/2.) + n) +
	   ((((1 + v.d - n)/(1 - v.m - n/2.) +
	           ((1 + v.d - n)*(1 + v.m - n/2.))/pow(1 - v.m - n/2.,2))/
	         (2.*sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.))) +
	        (-((v.d*(-v.m + n/2.))/pow(v.m + n/2.,2)) - v.d/(v.m + n/2.))/
	         (2.*sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.))))*
	      (sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.)) +
	        sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.)))*(-2*v.d + n))/(-2*(-v.m + n/2.)*(v.m + n/2.) + n);
}

// Using Mathematica
double Func::dgtdd() {
	return (-2*(sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.)) +
	        sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.)))*
	      (sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.)) +
	        sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.))))/(-2*(-v.m + n/2.)*(v.m + n/2.) + n) +
	   ((sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.)) +
	        sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.)))*
	      ((1 - v.m - n/2.)/
	         (2.*sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.))*(1 + v.m - n/2.)) +
	        (v.m + n/2.)/(2.*(-v.m + n/2.)*sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.))))*(-2*v.d + n))/
	    (-2*(-v.m + n/2.)*(v.m + n/2.) + n) +
	   (((1 + v.m - n/2.)/
	         (2.*(1 - v.m - n/2.)*sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.))) +
	        (-v.m + n/2.)/(2.*sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.))*(v.m + n/2.)))*
	      (sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.)) +
	        sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.)))*(-2*v.d + n))/(-2*(-v.m + n/2.)*(v.m + n/2.) + n);
}

double Func::mygamma(double kx, double ky) {
   return 2.*(cos(kx) + cos(ky));
}

double Func::eta(double kx, double ky) {
   return 2.*(cos(kx) - cos(ky));
}

double Func::epsilon(double kx, double ky) {
   return -(gt()*p.t + 3./8.*gs()*p.J*v.chi)*mygamma(kx,ky);
}

double Func::BDelta_d() {
   return 3./8.*gs()*p.J*v.BDelta;
}

double Func::BDelta_af() {
   return 2.*gs()*p.J*v.m;
}

/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	// result of f (independent on k_x and k_y)
	double gs;
	double dgtdm;
	double dgtdd;
	double dgsdm;
	double dgsdd;
	double BDelta_d;
	double BDelta_af;

	void onePoint(variables &y, int i, int j);

	void refreshRF() {
		gs = f->gs();
		dgtdm = f->dgtdm();
		dgtdd = f->dgtdd();
		dgsdm = f->dgsdm();
		dgsdd = f->dgsdd();
		BDelta_d = f->BDelta_d();
		BDelta_af = f->BDelta_af();
	}

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); refreshRF(); }
	~Equa() { delete f; }

	void setPar(parameters &p) {
		this->p = p;
		f->setPar(p);
		refreshRF();
	}

	void run(variables &w);
};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	double gamma = f->mygamma(x, y);
	double eta = f->eta(x, y);
	double epsilon = f->epsilon(x, y);
	double dzeta = sqrt(epsilon*epsilon + BDelta_af*BDelta_af);
	double a1 = dzeta-v.mu;
	double a2 = -dzeta-v.mu;
	double b = BDelta_d*eta;
	double E1 = sqrt(a1*a1 + b*b);
	double E2 = sqrt(a2*a2 + b*b);

	double part = (dzeta-v.mu)/E1 + (dzeta+v.mu)/E2;
	double dEvardt = p.t*gamma*epsilon/dzeta*part;
	double dEvards = (3./8.*v.chi*gamma*epsilon-2*v.m*BDelta_af)*p.J*part/dzeta -
			3./8.*p.J*v.BDelta*eta*eta*BDelta_d*(1./E1+1./E2);

	w.BDelta += eta*eta*BDelta_d*(1./E1 + 1./E2);
	w.chi += gamma*epsilon/dzeta*((-dzeta+v.mu)/E1 + (-dzeta-v.mu)/E2);
	w.m += BDelta_af/dzeta*part - (dEvardt*dgtdm + dEvards*dgsdm)/(2.*gs*p.J);
	w.d += dEvardt*dgtdd + dEvards*dgsdd;
	w.mu += (dzeta-v.mu)/E1 + (-dzeta-v.mu)/E2;

	w.Energy += E1 + E2;
}

void Equa::run(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation

	unsigned int max = p.L/2; // attention: L should be even!

	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
	// * edges (2 times)
	for(unsigned int i=1; i<max; i++) {
		onePoint(temp, i, 0);
		onePoint(temp, 0, i);
		onePoint(temp, i, max-i);
	}
	addTo(temp, w, 2);
	clear(temp);
	// * interior area (4 times)
	for(unsigned int i=1; i<max; i++) {
		for(unsigned int j=1; j<i; j++) {
            onePoint(temp, max-i, j);
		}
	}
	addTo(temp, w, 4);
	clear(temp);
	// * missing corners
	onePoint(w, max, 0);
	//onePoint(w, 0, max);
	onePoint(w, 0, 0);

	double chi2 = v.chi*v.chi;
	double m2 = v.m*v.m;
	double BDelta2 = v.BDelta*v.BDelta;

	double L2 = (double)(p.L*p.L);
	w.BDelta /= 4.*L2;
	w.chi /= 4.*L2;
	w.m /= 2.*L2;
	w.d /= L2;
	w.mu /= L2;
	w.Energy /= L2;

	w.BDelta += -v.BDelta;
	w.chi += -v.chi;
	w.m += -v.m - dgsdm*(3./4.*(BDelta2+chi2) + 2.*m2)/(4.*gs);
	w.d += p.U + dgsdd*(3./4.*(BDelta2+chi2) + 2.*m2)*p.J;
	w.mu += -p.sdelta;

	v.Energy = p.U*v.d - v.mu*p.sdelta - w.Energy + (3./4.)*gs*p.J*(BDelta2 + chi2) + 2.*gs*p.J*m2;

	//if(v.m < 0) w.m = v.m-0.4; // temporary for debuging
	if(w.BDelta != w.BDelta) { w.BDelta = 1000; }
	if(w.chi != w.chi) { w.chi = 1000; }
	if(w.m != w.m) { w.m = 1000; }
	if(w.d != w.d) { w.d = 1000; }
	if(w.mu != w.mu) { w.mu = 1000; }

	// results
	glob.Energy = v.Energy;
	glob.gt = f->gt();
	glob.gs = f->gs();
	glob.BDelta_sc = v.BDelta*f->gBDelta();
	glob.m_AF = v.m*f->gm();
	glob.BDelta_d = 3./8.*f->gs()*p.J*v.BDelta;
	glob.gBDelta = f->gBDelta();
	glob.gm = f->gm();
}

// ***********************************************************************************  //
// *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY *** END OF THE COPY ***
// ***********************************************************************************  //

int function(const gsl_vector * x, void *params, gsl_vector * f) {

	struct parameters p;
	p.sdelta = ((Par *) params)->sdelta();
	p.J = ((Par *) params)->J();
	p.t = ((Par *) params)->t();
	p.U = ((Par *) params)->U();
	p.L = ((Par *) params)->L();

	struct variables v;
	v.BDelta = gsl_vector_get (x, 0);
	v.chi = gsl_vector_get (x, 1);
	v.m =  gsl_vector_get (x, 2);
	v.d = gsl_vector_get (x, 3);
	v.mu = gsl_vector_get (x, 4);
	v.Energy = 0; // it's calculated later

	if(DEEP_DEBUG) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4));
	}

    // It's start of the calculations

    Equa E(p, v);
    variables w;
    E.run(w);

    // It's end of the calculations

    gsl_vector_set (f, 0, w.BDelta);
    gsl_vector_set (f, 1, w.chi);
    gsl_vector_set (f, 2, w.m);
    gsl_vector_set (f, 3, w.d);
    gsl_vector_set (f, 4, w.mu);

	if(DEEP_DEBUG) {
		printf("### %25.15e %25.15e %25.15e %25.15e %25.15e  E=%.16f\n",
				gsl_vector_get (x, 0), gsl_vector_get (x, 1), gsl_vector_get (x, 2),
				gsl_vector_get (x, 3), gsl_vector_get (x, 4), glob.Energy);
		printf("###  %24.15e# %24.15e# %24.15e# %24.15e# %24.15e#\n",
			gsl_vector_get (f, 0), gsl_vector_get (f, 1), gsl_vector_get (f, 2),
			gsl_vector_get (f, 3), gsl_vector_get (f, 4));
		printf("\n");
	}

    return GSL_SUCCESS;
}

// ************* CLASSES ************* //

// Par

void Par::print(FILE *file) const {
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.12e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
}

void Par::print_verboten(FILE *file) const {
	fprintf(file, "PARAMETERS:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.12e", tabD[i]); }
	for(int i=0; i<NtabI; i++) { fprintf(file, " %6d", tabI[i]); }
	fprintf(file, "\n");
}

void Par::print_label(FILE *file)  const {
	fprintf(file, "sdelta=%f, J=%f, t=%f, U=%f, L=%d", sdelta(), J(), t(), U(), L());
}

// Var

void Var::print(FILE *file) const{
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.12e", tabD[i]); }
}

void Var::print_verboten(FILE *file) const{
	fprintf(file, "VARIABLES:");
	for(int i=0; i<NtabD; i++) { fprintf(file, " %21.12e", tabD[i]); }
}

void Var::print_label(FILE *file) const {
	fprintf(file, "       BDelta                chi                   m                     d                     mu");
}

// Res

void Res::set(gsl_multiroot_fsolver *s) {
	for(int i=0; i<n; i++) {
		(*V)[i] = gsl_vector_get (s->x, i);
	}
	this->glob = case_1D1t_GA::glob;
}

void Res::print(FILE *file) const { V->print(file); fprintf(file, "%21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f %21.15f\n", glob.Energy, glob.gt, glob.gs, glob.gBDelta, glob.gm, glob.BDelta_sc, glob.m_AF, glob.BDelta_d); }
void Res::print_verboten(FILE *file) const { fprintf(file, "Result: "); V->print_verboten(file); fprintf(file, "E=%21.15f gt=%21.15f gs=%21.15f gBDelta=%21.15f gm=%21.15f BDelta_sc=%21.15f m_AF=%21.15f BDelta_d=%21.15f\n", glob.Energy, glob.gt, glob.gs, glob.gBDelta, glob.gm, glob.BDelta_sc, glob.m_AF, glob.BDelta_d); }
void Res::print_label(FILE *file) const { V->print_label(file); fprintf(file, "                      E                      gt                    gs             gBDelta               gm                    BDelta_sc             m_AF                  BDelta_d\n"); }


void example_inf() {
	Par par(N_PAR);
	par.set_sdelta(0.35);
	par.set_J(1./3.);
	par.set_t(1.);
	par.set_U(5.);
	par.set_L(1024);
	par.print_label(); printf("\n");

	// I'm choosing start point
	Var initV(N_VAR);

	// INFINITE
	Res res(N_VAR);
	RootGSL root(par, function, res);
	printf("sdelta      "); res.print_label();
	while(1) {
		initV.set_BDelta(frand());
		initV.set_chi(frand());
		initV.set_m(frand());
		initV.set_d(frand());
		initV.set_mu(frand()*16.-8.);
		if(root.Solve(initV) == GSL_SUCCESS) {
			printf("       "); root.print();
		}
	}
}

void example_1D() {
	Par par(N_PAR);
	par.set_sdelta(0.03);
	par.set_J(1./3.);
	par.set_t(1.);
	par.set_U(10.);
	par.set_L(256);
	printf("#"); par.print_label(); printf("\n");

	// 10.000000    -1.273626604528e-01    3.682367177240e-01    1.444275071372e-01    1.989872348413e-02   -1.513455412177e-01      -0.377328419446        0.248018053937        2.953448291073       -0.031786117645        0.248207328996       -0.047019878983
	// I'm choosing start point
	Var initV(N_VAR);
	initV.set_BDelta(1.273626604528e-01);
	initV.set_chi(3.682367177240e-01);
	initV.set_m(1.444275071372e-01);
	initV.set_d(1.989872348413e-02);
	initV.set_mu(-1.513455412177e-01);

	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);

	LOOPS = 100;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(10 + ((double)i)*0.1);
		printf("%f ", par.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 100;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(10 - ((double)i)*0.1);
		printf("%f ", par.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
}

void test_energy() {

	struct parameters p;
	p.sdelta = 0.6;
	p.J = 1./3.;
	p.t = 1.;
	p.U = 12.;
	p.L = 256;

//	struct variables v1;
//	v1.BDelta = 2.251021863958e-01;
//	v1.chi = 3.918123984194e-01;
//	v1.m = -4.265632387964e-14;
//	v1.mu = -3.794176833425e-02;
//	v1.d = 2.888477153509e-02;

	// faworyzowane
	// #sdelta  U              BDelta                chi                   m                     d                     mu                      E                      gt                    gs                    BDelta_sc             m_AF                  BDelta_d
	// 0.030000 10.000000    1.274575403202e-01    3.682224936904e-01    1.444495221848e-01    1.989267541031e-02   -1.512343887518e-01      -0.377360693582        0.247974013523        2.953386298534        0.031804296571        0.248242557843        0.047053919153
//	struct variables v2;
//	v2.BDelta = 1.274575403202e-01;
//	v2.chi = 3.682224936904e-01;
//	v2.m = 1.444495221848e-01;
//	v2.mu = -1.512343887518e-01;
//	v2.d = 1.989267541031e-02;

	// nierealizowane
	// #sdelta  U              BDelta                chi                   m                     d                     mu                      E                      gt                    gs                    BDelta_sc             m_AF                  BDelta_d
    // 0.030000 10.000000    8.844740952418e-02    3.649240650551e-01    2.533547287718e-01    2.985890790656e-02   -2.921101500515e-01      -0.365504866521        0.337597633706        2.101519589552        0.030128661283        0.367278694448        0.023234245470
//	struct variables v3;
//	v3.BDelta = 8.844740952418e-02;
//	v3.chi = 3.649240650551e-01;
//	v3.m = 2.533547287718e-01;
//	v3.mu = -2.921101500515e-01;
//	v3.d = 2.985890790656e-02;

    // It's start of the calculations
	/*
		GA
		sdelta             BDelta                chi                   m                     d                     mu                      E                      gt                    gs             gBDelta               gm                    BDelta_sc             m_AF                  BDelta_d
		           5.102236912890e-04    3.593609375214e-01    1.505490218941e-15    1.122392620654e-02   -9.029144426464e-01   -1.029980725110211     0.767764116249644     1.890966457949343     0.767764116249644     1.375124160921240     0.000391731441432     0.000000000000002     0.000120601985785

		SGA 1D1t_v1
		sdelta           BDelta                chi                   m                     mu                    d                     lBDelta               lchi                  lm                    ln                      E                      gt                    gs             gBDelta               gm
		          -6.924471251818e-13    3.593602178999e-01   -1.256790616034e-13    4.452148266829e-01    1.122388991608e-02    4.484767526434e-16   -4.795763648411e-15    3.828097687905e-13   -1.348310680017e+00   -1.029982378303645     0.767763802616342     1.890966933224481     0.767763802616342     1.375124333732947


		SGA 1D1t_v2
		sdelta          barchiAB             barDelta                mu                     mAF                    d                      E                      gt                    gs             gBDelta               gm                 lbarchiAB             lbarDelta             ln                  lmAF
		           3.593602179003e-01   -2.818628663568e-13    4.452148266826e-01   -3.646948119761e-13    1.122388991610e-02   -1.029982378303640     0.767763802616518     1.890966933224215     0.767763802616518     1.375124333732850     0.852706088762232    -0.000000000000067    -1.348310680017772     0.000000000001117
	*/
	struct variables v;
	v.BDelta = 6.924471251818e-13;
	v.chi = 3.593602178999e-01;
	v.m = 1.256790616034e-13;
	v.mu = 4.452148266829e-01 -1.348310680017e+00;
	v.d = 1.122388991608e-02;

    Equa E1(p, v);
    variables w1;
    E1.run(w1);
    printf("0_B> %25.15e\n", glob.Energy);

/*    Equa E2(p, v3);
    variables w2;
    E2.run(w2);
    printf("0_C> %25.15e\n", glob.Energy);

    Equa E3(p, v);
    variables w3;
    E3.run(w3);
    printf("0> %25.15e\n", glob.Energy);
*/
    // It's end of the calculations
    /*
    Func F(p, v);
    printf("1> %25.15e\n", F.BDelta_af());
    printf("1> %25.15e\n", F.BDelta_d());
    printf("1> %25.15e\n", F.dgsdd());
    printf("1> %25.15e\n", F.dgsdm());
    printf("1> %25.15e\n", F.dgtdd());
    printf("1> %25.15e\n", F.dgtdm());
    printf("1> %25.15e\n", F.gs());
    printf("1> %25.15e\n", F.gt());

    double x = (double)100*2.*M_PI/p.L;
    double y = (double)136*2.*M_PI/p.L;

    //double gamma = F.mygamma(x, y);
	double eta = F.eta(x, y);
	double epsilon = F.epsilon(x, y);
	double dzeta = sqrt(epsilon*epsilon + F.BDelta_af()*F.BDelta_af());
	double a1 = dzeta-v.mu;
	double a2 = -dzeta-v.mu;
	double b = F.BDelta_d()*eta;
	double E1 = sqrt(a1*a1 + b*b);
	double E2 = sqrt(a2*a2 + b*b);

//	printf("\n");
//	printf("1> x=%f, y=%f\n", x, y);
//    printf("1> E1     %25.15e\n", E1);
//    printf("1> E2     %25.15e\n", E2);
//    printf("1> E1+E2  %25.15e\n\n", E1+E2);
//    printf("1> Global %25.15e\t%f\n", glob.Energy, debug_e);

    printf("0> %25.15e\n", glob.Energy);
*/
    return;
}


void example_2D_U_run(Par &par, Var &initV, int s1, int s2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	int LOOPS_sdelta, LOOPS_U;
	Var initV_cp(N_VAR); initV_cp = initV;

	double start_U = 10;
	double start_sdelta = 0.01;

	par.set_U(start_U);
	par.set_sdelta(start_sdelta);

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta  U       "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	if(s2>0) { LOOPS_U = 400; } else { LOOPS_U = 200; }
	if(s1>0) { LOOPS_sdelta = 440; } else { LOOPS_sdelta = 10; }

	for(int i=0; i<=LOOPS_sdelta; i++) {
		if(s1>0) { par.set_sdelta(start_sdelta + (double)i*0.001); } else { par.set_sdelta(start_sdelta - (double)i*0.001); }
		par.set_U(start_U);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS_U; j++) {
			if(s2>0) { par.set_U(start_U + (double)j*0.05); } else { par.set_U(start_U - (double)j*0.05); }
			printf("%f %f", par.sdelta(), par.U());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
	}
	printf("\n#END\n");
	fflush(stdout);
	initV = initV_cp;
}

void example_2D_sdelta_run(Par &par, Var &initV, int s1, int s2) {
	Res res(N_VAR);
	RootGSL *root = new RootGSL(par, function, res);
	int LOOPS_sdelta, LOOPS_U;
	Var initV_cp(N_VAR); initV_cp = initV;

	double start_U = 10;
	double start_sdelta = 0.01;

	par.set_U(start_U);
	par.set_sdelta(start_sdelta);

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#U       sdelta        "); res.print_label();
	fflush(stdout);

	initV.shake();
	while(root->Solve(initV) != GSL_SUCCESS) { printf("# shot is not good.\n"); initV.shake(); return; }
	Var old(*res.get_V());

	if(s2>0) { LOOPS_sdelta = 440; } else { LOOPS_sdelta = 10; }
	if(s1>0) { LOOPS_U = 400; } else { LOOPS_U = 200; }

	for(int i=0; i<=LOOPS_U; i++) {
		if(s1>0) { par.set_U(start_U + (double)i*0.05); } else { par.set_U(start_U - (double)i*0.05); }
		par.set_sdelta(start_sdelta);
		initV = old;
		initV.shake();
		if(root->Solve(initV) != GSL_SUCCESS) { continue; }
		old.set(res.get_V());
		initV = old;

		for(int j=0; j<=LOOPS_sdelta; j++) {
			if(s2>0) { par.set_sdelta(start_sdelta + (double)j*0.001); } else { par.set_sdelta(start_sdelta - (double)j*0.001); }
			printf("%f %f", par.U(), par.sdelta());
			initV.shake();
			if(root->Solve(initV) != GSL_SUCCESS) { printf("\n"); fflush(stdout); continue; }
			root->print();
			fflush(stdout);
			initV.set(res.get_V());
		}

		printf("\n");
		fflush(stdout);
	}
	printf("\n#END\n");
	fflush(stdout);
	initV = initV_cp;
}

void example_2D() {
	Par par(N_PAR);
	par.set_sdelta(0.01);
	par.set_J(1./3.);
	par.set_t(1.);
	par.set_U(10.);
	par.set_L(512);
	printf("#"); par.print_label(); printf("\n"); fflush(stdout);

	// I'm choosing start point
	Var initV_A(N_VAR);
	//10.000000 0.03
	//initV_A.set_BDelta(2.251021892946e-01);
	//initV_A.set_chi(3.918123977692e-01);
	//initV_A.set_m(1.293183550440e-10);
	//initV_A.set_d(2.888477117053e-02);
	//initV_A.set_mu(-3.794176887261e-02);
	//10.000000 0.010000    2.441709755052e-01    3.877981484662e-01    4.918574899927e-13    2.915403194000e-02   -1.189892107432e-02   -0.283589086386898     0.253222268041979     3.472893999395024     0.253222268041979     1.863570229262912     0.061829528207453     0.000000000000917     0.105997489457311
	initV_A.set_BDelta(2.441709755052e-01);
	initV_A.set_chi(3.877981484662e-01);
	initV_A.set_m(4.918574899927e-13);
	initV_A.set_d(2.915403194000e-02);
	initV_A.set_mu(-1.189892107432e-02);

	Var initV_B(N_VAR);
	//10.000000 0.03
	//initV_B.set_BDelta(1.273626604528e-01);
	//initV_B.set_chi(3.682367177240e-01);
	//initV_B.set_m(1.444275071372e-01);
	//initV_B.set_d(1.989872348413e-02);
	//initV_B.set_mu(-1.513455412177e-01);
	//10.000000 0.010000    7.762016315035e-02    2.134028542536e-01    2.108432718951e-01    7.282345409992e-04   -1.709508969757e-01   -0.323201333655645     0.032677555467128     2.818174143172957     0.036763826270962     1.678741833389803     0.002853614193183     0.353951420819013     0.027343392097397
	initV_B.set_BDelta(7.762016315035e-02);
	initV_B.set_chi(2.134028542536e-01);
	initV_B.set_m(2.108432718951e-01);
	initV_B.set_d(7.282345409992e-04);
	initV_B.set_mu(-1.709508969757e-01);

	Var initV_C(N_VAR);
	//10.000000 0.03
	//initV_C.set_BDelta(8.844740599149e-02);
	//initV_C.set_chi(3.649240615837e-01);
	//initV_C.set_m(2.533547250583e-01);
	//initV_C.set_d(2.985890535207e-02);
	//initV_C.set_mu(-2.921101342172e-01);
	//10.000000 0.001000    1.776206018777e-02    3.490931987986e-01    2.800477162926e-01    3.103793586916e-02   -2.788334412306e-01   -0.249728093865128     0.262170978741626     2.034569554731782     0.262177865344469     1.426383382801336     0.004656819024149     0.399455408911285     0.004517268360918
	//initV_C.set_BDelta(1.776206018777e-02);
	//initV_C.set_chi(3.490931987986e-01);
	//initV_C.set_m(2.800477162926e-01);
	//initV_C.set_d(3.103793586916e-02);
	//initV_C.set_mu(-2.788334412306e-01);
	//10.000000 0.010000    5.484472731967e-02    3.525754693207e-01    2.742874005075e-01    2.896854847631e-02   -2.849065641094e-01   -0.289422209536476     0.276590657169023     2.053555617764506     0.277185059022450     1.433023243972165     0.015202138979173     0.393060220455913     0.014078337236509
	initV_C.set_BDelta(5.484472731967e-02);
	initV_C.set_chi(3.525754693207e-01);
	initV_C.set_m(2.742874005075e-01);
	initV_C.set_d(2.896854847631e-02);
	initV_C.set_mu(-2.849065641094e-01);

//	GA
//	sdelta=0.030000, J=0.333333, t=1.000000, U=10.000000, L=128
//	sdelta             BDelta                chi                   m                     d                     mu                      E                      gt                    gs                    BDelta_sc             m_AF                  BDelta_d
//	           2.251021892946e-01    3.918123977692e-01    1.293183550440e-10    2.888477117053e-02   -3.794176887261e-02      -0.368633521857        0.310899441561        3.334657313615        0.069984144946        0.000000000236        0.093829832730
//	           1.274575360262e-01    3.682224940937e-01   -1.444495202876e-01    1.989267564092e-02   -1.512343389056e-01      -0.377360693582        0.247974015166        2.953386307567        0.031804295702       -0.248242554962        0.047053917712
//	           8.844740599149e-02    3.649240615837e-01    2.533547250583e-01    2.985890535207e-02   -2.921101342172e-01      -0.365504866521        0.337597613555        2.101519638331        0.030128658318        0.367278693327        0.023234245081


/*	printf("# CASE A GA++\n"); example_2D_U_run(par, initV_A, 1, 1);
	printf("# CASE A GA+-\n"); example_2D_U_run(par, initV_A, 1, -1);
	printf("# CASE A GA-+\n"); example_2D_U_run(par, initV_A, -1, 1);
	printf("# CASE A GA--\n"); example_2D_U_run(par, initV_A, -1, -1);
*/
	//printf("# CASE B GA++\n"); example_2D_U_run(par, initV_B, 1, 1);
	printf("# CASE B GA+-\n"); example_2D_U_run(par, initV_B, 1, -1);
	//printf("# CASE B GA-+\n"); example_2D_U_run(par, initV_B, -1, 1);
	printf("# CASE B GA--\n"); example_2D_U_run(par, initV_B, -1, -1);
/*
	printf("# CASE C GA++\n"); example_2D_U_run(par, initV_C, 1, 1);
	printf("# CASE C GA+-\n"); example_2D_U_run(par, initV_C, 1, -1);
	printf("# CASE C GA-+\n"); example_2D_U_run(par, initV_C, -1, 1);
	printf("# CASE C GA--\n"); example_2D_U_run(par, initV_C, -1, -1);
*/

//	printf("# CASE A GA++\n"); example_2D_sdelta_run(par, initV_A, 1, 1);
//	printf("# CASE A GA+-\n"); example_2D_sdelta_run(par, initV_A, 1, -1);
//	printf("# CASE A GA-+\n"); example_2D_sdelta_run(par, initV_A, -1, 1);
//	printf("# CASE A GA--\n"); example_2D_sdelta_run(par, initV_A, -1, -1);

//	printf("# CASE B GA++\n"); example_2D_sdelta_run(par, initV_B, 1, 1);
//	printf("# CASE B GA+-\n"); example_2D_sdelta_run(par, initV_B, 1, -1);
//	printf("# CASE B GA-+\n"); example_2D_sdelta_run(par, initV_B, -1, 1);
//	printf("# CASE B GA--\n"); example_2D_sdelta_run(par, initV_B, -1, -1);

//	printf("# CASE C GA++\n"); example_2D_sdelta_run(par, initV_C, 1, 1);
//	printf("# CASE C GA+-\n"); example_2D_sdelta_run(par, initV_C, 1, -1);
//	printf("# CASE C GA-+\n"); example_2D_sdelta_run(par, initV_C, -1, 1);
//	printf("# CASE C GA--\n"); example_2D_sdelta_run(par, initV_C, -1, -1);
}

void run_1D_n() {
	Par par(N_PAR);
	par.set_sdelta(0.03);
	par.set_J(1./3.);
	par.set_t(1.);
	par.set_U(8);
	par.set_L(1600);
	printf("#"); par.print_label(); printf("\n");

		//#sdelta(0.01)    BDelta                chi                   m                     d                     mu                      E                      gt                    gs             gBDelta               gm                    BDelta_sc             m_AF                  BDelta_d
//SC	//12.000000     3.157343962587e-01    3.576283386987e-01   -2.731349105807e-18    2.689257842879e-03   -8.667579219219e-03   -0.264654660600024     0.053296556352707     3.878693569419671     0.053296556352707     1.969439912619745     0.016827556042690    -0.000000000000000     0.153079621551659
//SC+AF //12.000000     7.730712909198e-02    1.974978438726e-01    2.130523518337e-01    3.579267941554e-04   -1.726307983976e-01   -0.322198602105712     0.027256785313817     2.804530765330826     0.032179310506703     1.674673330930790     0.002487690111433     0.356793091708023     0.027101277739732
//SC+AF //12.000000     5.543852122622e-02    3.006964425218e-01    3.565405089760e-01    9.973337147538e-03   -3.071741396300e-01   -0.255210784623441     0.150007705232853     1.654348740833274     0.153567003306780     1.286214889057530     0.008513527572470     0.458587711197024     0.011464330973033

	// Najbardziej prawdopodobne
	//0.004000    -4.809055017588e-02    9.366616860808e-02    2.262798101456e-01    2.475691439082e-05   -1.799589427225e-01   -0.319755177506474     0.008276570848392     2.733441248770572     0.011153696545014     1.653312205474384    -0.000536387403345     0.374111171966119    -0.016431586690852

	// I'm choosing start point
	Var initV_A(N_VAR);
	initV_A.set_BDelta(3.157343962587e-01);
	initV_A.set_chi(3.877981484662e-01);
	initV_A.set_m(2.731349105807e-18);
	initV_A.set_d(2.689257842879e-03);
	initV_A.set_mu(-8.667579219219e-03);

	Var initV_B(N_VAR);
	initV_B.set_BDelta(-8.667579219219e-03);
	initV_B.set_chi(1.974978438726e-01);
	initV_B.set_m(2.130523518337e-01);
	initV_B.set_d(3.579267941554e-04);
	initV_B.set_mu(-1.726307983976e-01);

	Var initV_C(N_VAR);
	initV_C.set_BDelta(5.543852122622e-02);
	initV_C.set_chi(3.006964425218e-01);
	initV_C.set_m(3.565405089760e-01);
	initV_C.set_d(9.973337147538e-03);
	initV_C.set_mu(-3.071741396300e-01);

	Var initV_C2(N_VAR);
	initV_C2.set_BDelta(4.809055017588e-02);
	initV_C2.set_chi(9.366616860808e-02);
	initV_C2.set_m(2.262798101456e-01);
	initV_C2.set_d(2.475691439082e-05);
	initV_C2.set_mu(-1.799589427225e-01);

	//0.099000     1.940596348836e-01    3.941589456832e-01   -1.194162414833e-10    1.566740334787e-02   -1.295901349839e-01   -0.556602030314726     0.377843716342518     3.085454893026896     0.377843716342518     1.756546296864075     0.07
	Var initV_C3(N_VAR);
	initV_C3.set_BDelta(1.940596348836e-01);
	initV_C3.set_chi(3.941589456832e-01);
	initV_C3.set_m(1.194162414833e-10);
	initV_C3.set_d(1.566740334787e-02);
	initV_C3.set_mu(-1.295901349839e-01);

	// 0.03
	// 5.000000     7.806614184421e-02    4.044804316492e-01    1.103739334650e-14    1.298662678259e-01   -4.655831377911e-02   -0.765842647437563     0.821677624849884     2.021556648152390     0.821677624849884     1.421814561802062     0.064145202011741     0.000000000000016     0.019726891005095
	Var initV_B2(N_VAR);
	initV_B2.set_BDelta(7.806614184421e-02);
	initV_B2.set_chi(4.044804316492e-01);
	initV_B2.set_m(1.103739334650e-14);
	initV_B2.set_d(1.298662678259e-01);
	initV_B2.set_mu(-4.655831377911e-02);

	// duze U
	//SC+AF //12.000000     7.730712909198e-02    1.974978438726e-01    2.130523518337e-01    3.579267941554e-04   -1.726307983976e-01   -0.322198602105712     0.027256785313817     2.804530765330826     0.032179310506703     1.674673330930790     0.002487690111433     0.356793091708023     0.027101277739732
	Var initV_DuzeU(N_VAR);
	initV_DuzeU.set_BDelta(7.730712909198e-02);
	initV_DuzeU.set_chi(1.974978438726e-01);
	initV_DuzeU.set_m(2.130523518337e-01);
	initV_DuzeU.set_d(3.579267941554e-04);
	initV_DuzeU.set_mu(-1.726307983976e-01);

	//0.300000 8.000000    3.393285491580e-03    3.815175752385e-01   -2.183306353133e-14    2.919762580505e-02   -6.324909237552e-01   -1.031979713033357     0.781876067434048     1.988439334169601     0.781876067434048     1.410120326131639     0.002653128715838    -0.000000000000031     0.000843417792941
	Var initV_U8(N_VAR);
	initV_U8.set_BDelta(3.393285491580e-03);
	initV_U8.set_chi(3.815175752385e-01);
	initV_U8.set_m(2.183306353133e-14);
	initV_U8.set_d(2.919762580505e-02);
	initV_U8.set_mu(-6.324909237552e-01);

	// I'm choosing start point
	Var initV(N_VAR);
	initV = initV_U8;

	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);

/*	LOOPS = 30;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_sdelta(0.03 - ((double)i)*0.001);
		printf("%f ", par.sdelta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");
*/
	LOOPS = 210;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_sdelta(0.03 + ((double)i)*0.002);
		printf("%f ", par.sdelta());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
		fflush(stdout);
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
}

void run_1D_U() {
	Par par(N_PAR);
	par.set_sdelta(0.03);
	par.set_J(1./3.);
	par.set_t(1.);
	par.set_U(10.);
	par.set_L(256);
	printf("#"); par.print_label(); printf("\n");

	// I'm choosing start point
	Var initV_A(N_VAR);
	//10.000000 0.03
	//initV_A.set_BDelta(2.251021892946e-01);
	//initV_A.set_chi(3.918123977692e-01);
	//initV_A.set_m(1.293183550440e-10);
	//initV_A.set_d(2.888477117053e-02);
	//initV_A.set_mu(-3.794176887261e-02);
	//10.000000 0.010000    2.441709755052e-01    3.877981484662e-01    4.918574899927e-13    2.915403194000e-02   -1.189892107432e-02   -0.283589086386898     0.253222268041979     3.472893999395024     0.253222268041979     1.863570229262912     0.061829528207453     0.000000000000917     0.105997489457311
	initV_A.set_BDelta(2.441709755052e-01);
	initV_A.set_chi(3.877981484662e-01);
	initV_A.set_m(4.918574899927e-13);
	initV_A.set_d(2.915403194000e-02);
	initV_A.set_mu(-1.189892107432e-02);

	Var initV_B(N_VAR);
	//10.000000 0.03
	initV_B.set_BDelta(1.273626604528e-01);
	initV_B.set_chi(3.682367177240e-01);
	initV_B.set_m(1.444275071372e-01);
	initV_B.set_d(1.989872348413e-02);
	initV_B.set_mu(-1.513455412177e-01);
	//10.000000 0.010000    7.762016315035e-02    2.134028542536e-01    2.108432718951e-01    7.282345409992e-04   -1.709508969757e-01   -0.323201333655645     0.032677555467128     2.818174143172957     0.036763826270962     1.678741833389803     0.002853614193183     0.353951420819013     0.027343392097397
//	initV_B.set_BDelta(7.762016315035e-02);
//	initV_B.set_chi(2.134028542536e-01);
//	initV_B.set_m(2.108432718951e-01);
//	initV_B.set_d(7.282345409992e-04);
//	initV_B.set_mu(-1.709508969757e-01);

	Var initV_C(N_VAR);
	//10.000000 0.03
	//initV_C.set_BDelta(8.844740599149e-02);
	//initV_C.set_chi(3.649240615837e-01);
	//initV_C.set_m(2.533547250583e-01);
	//initV_C.set_d(2.985890535207e-02);
	//initV_C.set_mu(-2.921101342172e-01);
	//10.000000 0.001000    1.776206018777e-02    3.490931987986e-01    2.800477162926e-01    3.103793586916e-02   -2.788334412306e-01   -0.249728093865128     0.262170978741626     2.034569554731782     0.262177865344469     1.426383382801336     0.004656819024149     0.399455408911285     0.004517268360918
	//initV_C.set_BDelta(1.776206018777e-02);
	//initV_C.set_chi(3.490931987986e-01);
	//initV_C.set_m(2.800477162926e-01);
	//initV_C.set_d(3.103793586916e-02);
	//initV_C.set_mu(-2.788334412306e-01);
	//10.000000 0.010000    5.484472731967e-02    3.525754693207e-01    2.742874005075e-01    2.896854847631e-02   -2.849065641094e-01   -0.289422209536476     0.276590657169023     2.053555617764506     0.277185059022450     1.433023243972165     0.015202138979173     0.393060220455913     0.014078337236509
	initV_C.set_BDelta(5.484472731967e-02);
	initV_C.set_chi(3.525754693207e-01);
	initV_C.set_m(2.742874005075e-01);
	initV_C.set_d(2.896854847631e-02);
	initV_C.set_mu(-2.849065641094e-01);



	// I'm choosing start point
	Var initV(N_VAR);
	initV = initV_B;

	Res res(N_VAR);
	RootGSL root(par, function, res);
	int LOOPS;

	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);

	LOOPS = 100;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(10 + ((double)i)*0.1);
		printf("%f ", par.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	printf("#NOT END YET. I'm changing only direction. Still the same case\n\n\n");

	LOOPS = 100;
	initV = old;
	for(int i=0; i<=LOOPS; i++) {
		par.set_U(10 - ((double)i)*0.1);
		printf("%f ", par.U());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
}

void test_L() {
	Par par(N_PAR);
    par.set_sdelta(0.3);
    par.set_J(1./3.);
    par.set_t(1.);
    par.set_U(5);
    par.set_L(128);
	printf("#"); par.print_label();
	printf("\n#"); par.print_verboten();

	// #sdelta             BDelta                chi                   m                     d                     mu                      E                      gt                    gs             gBDelta               gm                    BDelta_sc             m_AF                  BDelta_d
	// 0.250000     3.118117003509e-03    3.895642522910e-01   -7.316060618842e-15    6.187739037283e-02   -5.436728568159e-01   -1.114826437637910     0.870518140428657     1.784868748578874     0.870518140428658     1.335989801075919     0.002714377415534    -0.000000000000010     0.000695678699247
	// 0.300000    -5.726209623548e-04    3.815190662602e-01   -6.750517100808e-17    5.222685987114e-02   -6.969481175397e-01   -1.150929046162602     0.884438893719889     1.713200685563590     0.884438893719889     1.308892923643332    -0.000506448250466    -0.000000000000000    -0.000122626828159
	Var initV_B2(N_VAR);
	initV_B2.set_BDelta(3.118117003509e-03);
	initV_B2.set_chi(3.895642522910e-01);
	initV_B2.set_m(1.103739334650e-14);
	initV_B2.set_d(6.187739037283e-02);
	initV_B2.set_mu(-5.436728568159e-01);

	Var initV = initV_B2;

	Res res(N_VAR);
	RootGSL root(par, function, res);

	initV.shake();
	if(root.Solve(initV) != GSL_SUCCESS) { printf("First shot is not good.\n"); return; }
	Var old(*res.get_V());

	printf("#init_p   "); initV.print(); printf("\n");
	printf("#sdelta      "); res.print_label();
	fflush(stdout);


	int LOOPS1 = 10;
	int LOOPS2 = 10;
	int LOOPS3 = 10;
	int LOOPS4 = 10;
	int LOOPS5 = 10;
	int LOOPS6 = 10;

	initV = old;
	for(int i=0; i<LOOPS1; i++) {
		par.set_L(128 + 16*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<LOOPS2; i++) {
		par.set_L(288 + 32*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<LOOPS3; i++) {
		par.set_L(608 + 64*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<LOOPS4; i++) {
		par.set_L(1248 + 128*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<LOOPS5; i++) {
		par.set_L(2528 + 256*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}
	for(int i=0; i<=LOOPS6; i++) {
		par.set_L(5088 + 512*i);
		printf("%d ", par.L());
		initV.shake();
		if(root.Solve(initV) != GSL_SUCCESS) { printf("\n"); initV = old; fflush(stdout); continue; }
		root.print();
		initV.set(res.get_V());
	}

	printf("#END OF THE CASE\n\n\n");
	fflush(stdout);
}

} // end of namespace case_1D1t_GA
